February 1, 2008 
Tashkent 



Complex Diffusion Monte-Carlo method for the systems with 

complex wave function: 
test by the simulation of 2D electron in uniform magnetic field 



B. Abdullaev 1 , M. Musakhanov 1 , A. Nakamura 2 
1 Theoretical Physics Dept, Tashkent State University, 
Tashkent 700174, Uzbekistan 
e-mail: bah_abd@iaph.silk.org, yousuf@iaph.silk.org 
2 RUSE, Hiroshima University, Japan 



Abstract 

On the base of Diffusion Monte-Carlo method it is developed a new Complex 
Diffusion Monte-Carlo (CDMC) method allowing to simulate the quantum 
systems with complex wave function. There are no approximations on the 
calculation of modulus and phase of wave function in contrast to other meth- 
ods. We find that the averaged value of any quantity in CDMC will have 
no direct contribution from the phase of the distribution function but only 
from the phase of the Green function of the diffusion equation. This is most 
important and crucial point of CDMC. 

We are testing CDMC by the calculation of the wave function and the 
ground state energy of two-dimensional electron placed into the external uni- 
form magnetic field. There is an excellent agreement between simulations' 
results and an analytical ones. 

I. INTRODUCTION 



Quantum-mechanical wave function is an essentially complex in many cases. There are 
well-known examples - electron in the external uniform magnetic field when vector potential 
has a central-symmetric form and system of anyons. The complexity of wave function does 
not give a possibility for the simulations of such kind of systems by using well-known Green 
Function Monte-Carlo method (see review [p]]). This method essentially demand the reality 
of the system's wave function which considered as probability weight in during stochastic 
process. 

There had been undertaken the several attempts to construct a Monte-Carlo method 
for the simulation of quantum systems with complex wave function. Authors of the works 
suggested to take as a probability weight the module of the complex distribution func- 
tion. All quantities are calculated in |2||| by averaging over this complex distribution func- 
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tion. The phase of this function a is taken accordingly a = J2 a i, where aj is the phase at 
the z-time step of the stochastic process. 

Another quantum Monte Carlo method using algorithm without branching for the sim- 
ulating of complex problems was developed in [|J. The main trouble of this approach was 
an increasing of the statistical error as a function of the whole time of the simulation. 

The basic difficulty of the numerical simulation of the fermions is essentially the same as 
for the systems with the complex wave function. Their wave function can change the sign 
and therefore can not be used as the probability weight in the simulation process. 

For the simulation of the continuum (not on the lattice) fermionic systems it were de- 
veloped widely used fixed node Monte-Carlo method ||( see also reviews and recently 
proposed constrained path Monte-Carlo method |§. In these both methods it was assumed 
the restriction on the random walks connected with the uncertainties in the space local- 
ization of the wave function node surfaces. The comparison of these methods was done in 

E- 

Our method is very close to fixed phase Diffusion Monte-Carlo method developed in ||, 
which was applied to two-dimensional electrons in magnetic field. Laughlin's wave function 
T0|| was used as a trial one. In the framework of Diffusion Monte-Carlo (DMC) method 



it was calculated only modulus of the system's wave function. The phase of system's wave 
function was not calculated but considered as fixed and equal to the phase of Laughlin wave 
function. 

Here it is proposed new Complex Diffusion Monte-Carlo (CDMC) method including also 
the simulation of the phase factor of the wave function. 

CDMC are tested first by the calculations of the ground state wave function and energy 
of two-dimensional electron placed into the external uniform magnetic field. There is an 
excellent agreement between simulations' results and an analytical ones. Namely, it is re- 
produced the ground state energy and his degeneracy on the orbital quantum number m and 
also the simulated wave function is in exact correspondence with analytical prediction. The 
description of CDMC method is given at the section II of the paper. In the section III we 
represent the description of the simulation algorithm. Analytical expressions for the algo- 
rithm quantities are given at section IV. Section V represent the discussion of the simulation 
results. 



II. COMPLEX DIFFUSION MONTE-CARLO METHOD FOR THE QUANTUM 
SYSTEMS WITH COMPLEX WAVE FUNCTION 

Hamiltonian of the electron placed into external uniform magnetic field with vector 
potential A = - [Hr\, where H - magnetic field, r - electron's radius vector, and in external 
potential V(r) has a form: 

H = ^{p+^Af + V{r). (1) 

Here M - mass and e = — e charge of electron, p = —ihV, where V = i— — h j-^—, c - light 

ox Oy 

velocity. By taking in to account a relation pA — Ap = divA = for the such kind of vector 
potential one can rewrite Hamiltonian ([!]) in the following form: 
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Let us introduce a complex distribution function 

f(r,t) = ** T (f f )*(r,t). (3) 

Here ^(r) is a complex conjugated trial wave function of electron in magnetic field. Wave 
function t) satisfies a Schrodinger equation with imaginary time ( expressed in % units 

) 

-^§^ = (H-E T )V(r,t). (4) 

It is necessary to choose a trial energy in such manner || that at t — > oo \l/(r, £) — > 
\I/o(r) " an exac t stationary wave function of ground state of Hamiltonian (jl]). 

So, by accounting (f|), one can write the equation for the distribution function / = /(r, £): 

' ;y /7 V - + + c^(*o - ^)/, (5) 



dt 2M J M w vv " Mc 



where 



F q {t) = ^ T -\r)^* T (f), (6) 
E L (f) = y* T - l (r)H'y* T (r), (7) 

n 2 

For the p = —ihV and D = , and by renaming 

2M 

F Q (r) = 2** T - 1 (r)V** T (r), (9) 

we have an equation: 

" ^ = ~° Af + D ^ fP> *W) ~ ^ A ^f + - ( 10 ) 

Here A = V 2 . 

In general case when ^(r) is a complex wave function the quantity Fq(r) has a form: 

F Q (f) = ReF Q (f) + HmFQ (f) . (11) 
So, a final view of an equation for the distribution function / is 

= —DAf + DV(f ReF Q (r)) + i[V(Df ImF Q {f))- 



h\e\ 
Mc 



AVf] + (E L (r) - Er)f. 



Following 0, we assume: when the time step of integration of equation (|T2|) r — > the 
function Fq(r) remains constant, i.e. we will suppose that Fq(f) = Fq(f f ) at r — - > 0, where 
vector r corresponds to time point t + r and r ' to point t. 

Let us introduce a new quantity 



An(f, f') = -ImFnir ') — A(r). 

yv ' ' 2 QK ' 2DMc K ' 



(13) 



Then at r — > the Green function of equation ([12]) has a form: 
G(r, r '; r) = Gi(r, r r) exp [— t{El{t) — Et)\ exp iAg(r, f '){r — f ' — DrReFqif ')) 



where 



Gi (r,r ; t) = ^_ ^_ exp 



47rDr 



' - f ' - DrReF Q (r ')f 



(14) 



(15) 



One can see that the Green function G(r,f'; r) ( |i~4"D and distribution function / given by 
(|) are the complex functions. These both quantities are related by usual integral equation: 



/(r, t + r)= Jdf >G(r, f '; r)/(r i). 



(16) 



From equation (|lq) it is followed that the modulus and the phase of the distribution function 
at consequent time point are determined by the modulus and the phase of the Green function 
and of ones of the distribution function at previous time point of integration of diffusion 
equation ([T2|) . 

We have to note that for the general case of complex wave function ^(r*) the energy 
Eiif) is also a complex, so real and imaginary part of Ei(f) contribute to ones of the Green 
function. Therefore the last two exponents in ( |I4T) has a form: 



x exp 



exp [-r(ReE L (r) - E T )\ x 
iAq(r, f ')(r — f ' — DrReFqif ')) — irImEL(r) 



(17) 



III. THE DESCRIPTION OF THE SIMULATION ALGORITHM 

1) Let us to have the N c of the initial configurations r, i.e. a set of initial systems 
in which the electron position has random and uniform distribution. One can choose the 
configurations of r and with the distribution function /(r, 0) = |\I/t(^)| 2 , because at t = 
(0) is real. 

The choice of the boundary condition depends on the problem. For the one electron in 
the external magnetic field a boundary is periodic, i.e., if electron cross a boundary from 
one side of simulated cell, it enters into cell from opposite side. For example, for the system 
of bosons in 2D parabolic well and in external magnetic field the boundaries must be free 
because parabolic well determines itself the boundary of space distribution of the particles. 
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2) It is performed the quantum drift and diffusion of the particle from k - th configuration, 
for example, in accordance with the formula 

f k = f k ' + DrReF Q (f k ') + X . (18) 

Here x ~ 1S a gaussian random number having mean value zero and dispersion 1\J Dr. 

3) The transition into new space point in this configuration is accepted with probability 

P(f ' ->f,r) =min(l, W(r,f ')), (19) 

where 

W(f,f') 



|*t(0 | 2 Ci (r',r,r) 



|^T(r ')| 2 Gi(f,r ';r)' 

Here Gi(r, f r) is given by (|TJ). 

If transition of electron is accepted then in accordance with ( |17D it has new phase, if no 
then electron keep his old phase. 

4) After changing of the electron position from k - configuration into new space point, 
are calculated ReEi(r k ), ImE^ifk) and other quantities of interest. 

5) By using of the first exponential factor in (|lj]) it is calculated the multiplicity M k (the 
branching probability) for the configuration k accordingly 

M k = exp[-r(ReE L (f k ) - E T )\ (20) 

If M k is not integer, we add an uniformly distributed random number between and 1 to it 
and take M k equal to nearest integer. 

6) If M k ^ then M k copies of new fc-th configuration place in the list of the new 
N configurations that one is the initial at the next step r of integration of the diffusion 
equation. If M k = then there is no k-th configuration in the list of new iV configuration 
for the next time step r. 

All quantities of interest as ReEL(f k ), ImE^if^) and so on are multiplied by the factor 
M k for the calculating of mean values of these ones. 

7) It is repeated steps 2)-6) of algorithm until all N c configurations will not overlooked 
and electrons on these configurations will not simulated on the displacement and the having 
a new phase. 

8) It are calculated a mean energy and others mean quantities on the N number of 
configurations, got in the point 6) of algorithm, at this time step r in accordance with formula 
(p3|) (see below, where it is necessary to change M by N and to take into account that the 
phase of electron is equal to the phase of the Green function «g and at the calculation of 
a mean energy that the energy EL(r k ) has a real ReEL(r k ) and an imaginary ImE^rk) 
parts.) 

9) It are repeated steps l)-8) an integer number of time steps r of integration of diffusion 
equation. After that it are determined the mean values of quantities ReE, ImE and others 
on this integer number too. It is redetermined the new value of Et in accordance with 
(E^new = [(ET)oid + ReE]/2 accordingly with assumption that ImE « ReE. An integer 
number of time steps r represents an one time block At. 

10) Every time block At has the N c of initial configurations. The iV c of initial configu- 
rations are filled randomly by configurations of just ended time block. The random choice 
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of configurations consists of two steps: 

a) a random choice of a number of time step r in every time block At; 

b) a random choice of configuration from iV configurations at this fixed time step r. 

In this manner filled list of N c configurations will be initial one for the next time block At. 

11) The repeating of big number of time blocks At decreases essentially a correlation 
between configurations in neighbor time blocks and provides a right calculation of mean 
quantities. 

Let us discuss in detail the calculation of the mean quantities. In general case (see [@,0) 
the mean value of the some quantity F(R(t)) with the complex distribution function f(R, t) 
at the time t of the running process is 

M 

£exp[ia(i^(t)]F(%)) 

< F(t) >= . (21) 

^2exp[ia(Ri(t)] 

Here Ri{t) is the coordinates ri,r2, ...,fjv, N - number of particles of system, M - number 
of configurations at time point t. The particles coordinates fi, r*2, fjv in (21) are weighted 
with probability \f(R i ,t)\, and the quantity a(Ri(t)) is a phase of the distribution function 

f(Ri,t)- 

As it is clear from the integral expression flT6|), the phase of the distribution function at 
consequent time moment determines through the phase of the Green function and the one 
of the distribution function at previous time moment. So, we have: 

a(Ri(t + r)) = a G {Ri{t + r), R^t)) + a(Rj(t)). (22) 

Here ac(Ri(t + r), Rj(t)) is the phase of the Green function (fH]) (with account of (|17D), and 
index j shows that Rj taken from configurations at time moment t. By substituting (22) 
into fl2"ip we have 

M 

y^ e i QG (fl i (t+r),iij(t))+ia(iij(t)) F(R.(t _|_ r )) 

< F{t + r) >= i=i 



M 

g^G(« l (t+T),-R J (t))+ia(/? J (t)) 
i=l 

M 

e ia G (Ri (t+r),Rj (*)) _p ^ ^ + r )) 

4 = 1 

M 



(23) 



E 

8=1 

The mean quantity < F(t+r) > in (^) is determined only by the phase a G of the Green func- 
tion, because all quantities under sum in the numerator and the denominator are weighted 
with the probability \f(Rj,(t +t), t + t)\, i.e. at time moment t + r with new configurations 
Ri(t + r). 

Next, let us consider the expression for the phase of the Green function (0). We have 



6 



ac = Aq(f,f ')(r — r' — DrReFq{r ')) — rImEi{r). 



(24) 



From (|) and (0) it is seen that the expressions ReFq(f')) and El(t) have no r dependence. 
From (|18|) one can show that at r — ► \ f\ — ► |r '| as r 1 / 2 , because x has r 1 / 2 dependence. 
At the same limit Aq(r,r') — > Aq(f',r ') (see the expression ( |i~3|) for the Ag(r, r ')), i.e. 
there is no r dependence, also. So, we have «g ~ t 1//2 - 



IV. ANALYTICAL EXPRESSIONS FOR THE ALGORITHM QUANTITIES 

We take trial wave function of electron in the form: 

C f ax + t(3y \ m ( (x 2 + y 2 ) \ 

$j (r) = — exp -7— . (25) 

a H \ a H J \ Aa H ) 

Here C is normalization constant; a# = {%/ Mwh) 1 ^ 2 is magnetic length, where wh = 
\e\H/Mc - cyclotron frequency; x and y the two spatial components of electron's coordinate; 
m is orbital quantum number; a, (5 and 7 are arbitrary numerical constants. 

Ata = /3 = 7 = 1 (|25|) gives an exact analytical expression of electron ground state wave 
function. We deform exact electron wave function, i.e. take as the input a, j3 and 7 not 
equal to 1, and then study a relaxation of the trial function to the exact one via simulation 
process. 

As length unit we take an, as energy unit - hwn/2, and as time unit - 2/Hwh- 
By substituting of the wave function (^) into (H), we have 

ReF Q (r) = 1 [i(2ma 2 x - ja 2 x 3 - ^0 2 xy 2 ) + 
a z x z + p z y z 

+](2mP 2 y - ja 2 yx 2 - 7/?V )], (26) 

ft 2maf3 
ImF Q {r) = [-yi + x 3 . 

a z x z + p z y z 

Here and below i and j the unit vectors in x and y directions correspondingly. 
For the Aq(r,r') (|T3|) gives 



Mr, r ') = aV ^ v ^ (-2ma/V + y(aV 2 + p 2 y' 2 ))+ ^ 
+j(2maf3x' - x(a 2 x' 2 + f3 2 y' 2 ))}. 

The wave function fl25|) and eq. (0) (here V(r) = 0) gives 

#e£ L (r) = ReE^ (f) + ReE^ (f) , (28) 



where 



ReE^\f) = 7 (m + 1) + (x 2 + y 2 ) 



1 — 7 2 ma/3 



4 ra 2 + /3 2 y 2 



faEflr) = m ( m ~ ^(z 92 ~ « 2 )(^ 2 ~ gV) , 



(a 2 x 2 + (3 2 y 2 ) 2 
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and also 



ImEAr) 



mxy(j3 2 — a 2 ] 
a 2 x 2 + j3 2 y 2 



(m - l)2af3 
a 2 x 2 + j3 2 y 2 



(29) 



By using formulas (|2 
tities. 



9D , it is easy to write the analytical expressions for other quan- 



V. THE DISCUSSION OF THE SIMULATION RESULTS. 

The simulation of the ground state of one electron in magnetic field is a simple problem 
and was performed on the usual personal computer. 

Figs. 1-4 are presented the result of the calculations of the real ReE and imaginary part 
ImE of the ground state energy for different orbital quantum numbers m and different 
starting trial functions defined by the parameters a, j3, 7 as a functions of the number T of 
time blocks At (see point 9) of algorithm) ( a full number of time blocks are chosen ten in 
each running). 

In every time block At initial number of configurations N c was chosen equal 1000 and 
number of time steps r equal 100. We everywhere suggested r equal 0.01. Typical mean 
number of the population number N ( see point 6) of the algorithm) was around 2500 at 
initial blocks At and around 1000 at final blocks At. 

Fig.l present ReE and Fig.2 present lOOOImE for m = 0, 4, 8 with parameters a = (3 = 
7 = 1 as a functions of T. 

It is seen from Figs. 1,2 ReE is approaching fast to the exact value of ground state energy 
and reproduce very well an orbital quantum number m degeneracy, while the imaginary part 
of energy ImE ~ 10~ 3 ReE. 

The additional calculations shows, that the increasing of the statistics, i.e. the increasing 
of N c and a number of steps r in each block At, decreases ImE more. 

The Figs. 3, 4 represent the ReE and ImE for the different choices of the m = trial 
function. From these Figs, we can see that the deformation of wave function in some limits 
does not alter of the relaxation of the running process to exact final state. 

Fig. 5 present the initial uniform spatial distribution of electron before simulation, while 
Fig. 6 - the final spatial distribution of electron after simulation for m = 13. 

In this simulation the trial wave function ^ 1S taken with parameters a — [3 — 7 = 1, 
m — 13 and N c = 500. Since the trial function coincide with exact wave function we have 
to expect fast convergence of initial distribution to exact one if the algorithm is correct. In 
accordance with analytical solution (see [|TTJ) the spatial distribution for one electron in this 
case must have a ring like form with mean radius (2m + 1) 1 / 2 and wide 1 ( in magnetic length 
as units). As it is seen from Fig.6, there is an excellent agreement between simulation and 
analytical results. 



We suppose also further tests of the method by simulations of anyons [12 . 
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FIG. 1. The real part of the ground state energy ReE as a functions of the number T of time 
blocks At (see point 9) of algorithm) for different orbital quantum numbers m. Solid line - m = 0, 
marked solid line - m = 4, dashed line - m = 8. 
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FIG. 2. The same as Fig.l but for WOOImE. 
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FIG. 3. The real part of the ground state energy ReE for different starting parameters a, (3, 7 
and the same m = as a functions of the number T of time blocks At. Solid line -a = /3 = 7 = l, 
marked solid line - a = 1.4, (3 = 7 = 1, dashed line - a = /3=l,7 = 0.9. 

4. 00 q 

3. 00 ^ 

2. 00t 
UJ = 
£ j 
1—1 1 . 00- 

s> I 

q 0.00- 

- 1 . 00t 
-2. 00 

-3.00 1 

0.00 2.00 A. 00 6.30 8.00 10.00 12.00 

T 




FIG. 4. The same as Fig.3 but for WOOImE. 
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FIG. 5. Initial spatial distribution of the electron. 
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FIG. 6. The result of the simulation: final spatial distribution of the electron in uniform mag- 
netic field at m = 13. The length unity is an, initial number of points N c = 500. 
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